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Abstract 


The process of diffusion is the most elementary stochastic transport process. Brow¬ 
nian motion, the representative model of diffusion, played a important role in the ad¬ 
vancement of scientific fields such as physics, chemistry, biology and finance. However, 
in recent decades, non-diffusive transport processes with non-Brownian statistics were 
observed experimentally in a multitude of scientific fields. Examples include human 
travel [l[[2], in-cell dynamics [3], the motion of bright points on the solar surface [I], 
the transport of charge carriers in amorphous semiconductors [5], the propagation of 
contaminants in groundwater |6 , the search patterns of foraging animals [7 and the 
transport of energetic particles in turbulent plasmas [8 13 . These examples showed 
that the assumptions of the classical diffusion paradigm 14-16, assuming an un¬ 


derlying uncorrelated (Markovian), Gaussian stochastic process, need to be relaxed to 
describe transport processes exhibiting a non-local character and exhibiting long-range 
correlations. 

This article does not aim at presenting a complete review of non-diffusive trans¬ 
port, but rather an introduction for readers not familiar with the topic. For more in 
depth reviews, we recommend the references 17-20 . First, we recall the basics of 


the classical diffusion model and then we present two approaches of possible general¬ 
izations of this model: the continuous time random walk (CTRW) and the fractional 
Levy motion (fLm). 
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1 Classical diffusion and the random walk model 


The model of the random walk was first developed by Einstein [14| in 1905 and, inde¬ 
pendently, by Smoluchowski 115] in 1906, to explain the observation made by a Scottish 
botanist, Robert Brown (1773-1858), of the random motion of pollen particles in sus¬ 
pension on water. The term Brownian motion has been coined in honors of Brown to 
designate the random walk. In 1913, in his book Les Atomes |21|, Perrin verified the 
results of Einstein and Smoluchowski by measuring with a microscope the displacement 
of small particles in suspension in a liquid. Figure [l] reproduces some of his observations. 
Using Einstein’s theory, he was able to measure the Avogadro Number 22 . This remark¬ 


able success bore the definite proof of the existence of the atom which awarded Perrin the 
Nobel Prize for physics in 1926. 


Einstein’s model of the random walk assumes that each individual particle motion is 
independent of the other’s and that the displacements of the same particle at different 
times are also independent, provided that the interval separating the different times is not 
too small. Next, he introduces a time interval, r, very small compared to the observation 
time but sufficiently large so that the motions between two consecutive time intervals r can 
be considered as independent. In a one-dimensional (ID) model, considering n particles, 
during the time interval r, each particle position along the x-axis will increase by a value 
A, different for each particle. The probability distribution function (PDF) of the step 
sizes (p( A) is defined by the following relation: the number dn of particle experiencing a 
displacement lying between A and A + d A is given by 


dn = A) d A. 


(1) 


The step size probability satisfies the relation y?(A) = ip(— A) and differs from zero only 
for very small values of A. 


Let f(x,t ) be the distribution of particles. Using eq. [TJ the number of particles at time 
t T r found between x and x + dx is written 


f(x, t + t ) dx 


f + OO 


dx 


f(x + A, t)<p( A) d A. 


(2) 


Since r is very small, we can write 


9 f 

f(x,t + r) = f(x,t)+r - 7 ^, 
and by expanding f(x + A, t) in powers of A, we find 


. . . 9 f(x,t) A 2 9 2 f(x,t) . 3 . 

f(x + A, t) = /(*, t) + A + 0( A 3 ) 


d 2 x 


( 3 ) 

( 4 ) 


We note that, after the assumption of independent time steps and symmetrically dis¬ 
tributed step sizes, we assume here that r and A cannot take large values in order to 
perform the two previous expansions. This is, as we will see below, a fundamental restric- 
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Figure 1 | Example of trajectories of small particles (radius of 0.53 pm) in suspension in a fluid 
measured by Jean Perrin 23 . The successive positions are marked every 30 s and joined by straight 


lines to guide the eye. The grid size is 3.2 pm. 


tion of the diffusive model which gives it its local character in time and space. 


Since only very small values contribute to it, the expansion can be performed under the 
integral. We find 

Bf r+°° af f+°° B 2 f r +oc A 2 

fM+rjZ-fJ ,(A) d A + ^y Aip(A) dA+-^~ J — <p(A) dA+0 (A 3 ). 

(5) 

All the terms with odd powers of A vanish due to the fact that cp( A) = </?(—A). Taking 
into account the fact that J'+°° <p(A) d A = 1 , defining the variance of the step sizes 


/ +oo 

AV(A)dA 

-OO 

and retaining the terms up to 0(A 3 ), we find the well-known diffusion equation 

9/ = d 9V 

9 t 9 2 x’ 


( 6 ) 


(7) 


where D = ' 2 ' is the diffusion coefficient. The equation of diffusion had already been 
discovered experimentally by Fick in 1855, but Einstein was the first to derive it from this 
physical and atomistic model. 


Let us now find the fundamental solution, G(x, t), of eq. [7J i.e. the solution of the equation 
with initial condition G(x,t = 0) = 5(x). This solution is also called the Green function 
or the propagator of the equation. Its convolution with an arbitrary initial condition fo(x) 
provides the evolution of the initial condition at all times t > 0 

/ +oo 

G(x - x',t) fo(x') dx'. (8) 

-OO 
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As we will see later, it is interesting to take the Fourier transform in space of eq. [7] 

9 G(k,t) 


9 t 


= —kDG(k, t), 


where 


/ +oo 

-OO 


e~ ikx G(x,t) dx. 


The solution of eq. [9] i 


is 


G(k,t) = Ce 


—k 2 Dt 


( 9 ) 

( 10 ) 

( 11 ) 


As G(k, 0) = f_™ e lkx 5(x) dx = 1, we find that C = 1. 

We find the solution by taking the inverse Fourier transform of the previous expression 




p+oo -| /*+oo 

e - k 2 Dt+ikxdk = _ 

-no 27T 


— [ e d k. (12) 

After a little effort, one finds that the solution is the Gaussian, or normal , distribution 

(13) 


1 X 2 

G(x,t) = —j==.e iot, 


\J AnDt 

with a variance, or mean-squared displacement, given by 

(x 2 ) = 2 Dt. 


(14) 


The linear time dependence of the mean-square displacement of diffusive processes is a 
fundamental results of the random walk model. 

The fact that the distribution of positions of the random walkers is a Gaussian distribution 
arises naturally from the central limit theorem (CLT). Indeed, the position of each particle 
is a sum of independent and identically distributed steps having the same mean and the 
same variance. Each step being independent, the diffusive process is also memory-less, i.e. 
Markovian. 


2 Non-diffusive transport 


Since its discovery, the model of diffusion and Brownian motion plays a crucial role not 
only in physics but also in biology, chemistry, sociology, economics and finance. How¬ 
ever, numerous examples of transport show a deviation from the diffusive paradigm given 
by eq. (14). Non-diffusive transport is in fact ubiquitous in nature. Examples ranges 
from the dispersal of bank notes |1|, the motion of particles inside living cells |3| or the 
foraging movements of spider monkeys [7]. In plasma physics, examples are also multi¬ 
ple. The acceleration of electrons and ions by interplanetary shocks in the solar wind [8], 
ion transport across the magnetopause 24 , the motion of magnetic bright points on the 
solar surface j4| or the transport of tracer particles in 3D pressure-gradient driven turbu- 
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lence 125 are evidences indicating the presence of non-diffusive transport in magnetized 
plasmas. 


Non-diffusive transport is characterized by a mean-squared displacement (variance of dis¬ 
placement) of an ensemble of individuals that does not necessarily scale linearly with time 


(r(t) - {r(t ))) 2 ^} oc f 7 , 


(15) 


where r(t ) represents the positions of individuals and 7 is called the transport exponent. 
When 7 > 1 or 7 < 1, the transport is called superdiffusive or subdiffusive, respectively. 
For the special case of classical diffusion, 7 = 1 in accordance with eq. (14). When 


7 = 2, the transport is ballistic. Non-diffusive transport is at the heart of many complex 
systems, such as turbulence, where well-defined scale-lengths or time-scales do not exist, 
and thus transport cannot be modeled as a classical diffusive process. In these systems, 
the transport is characterized by the presence of long-term memory and/or non-Gaussian 
(heavy-tailed) PDFs. Mandelbrot coined the terms Noah effect and Joseph effect 26 , as 


a reference to the natural events experienced by these biblical figures, to describe those 
two effects in the context of hydrology. The great flood experienced by Noah and the 
seven years of abundance followed by seven years of famine experienced by Joseph are well 
known examples that reflects that extreme events with low probability and cycles or trend 
do, in fact, occur in nature. 


In order to account for these effects, the hypotheses of the CLT need to be loosen. By 
removing the restriction on the finiteness of the variance of the random variables (here, 
the step sizes), we allow large fluctuations in the random walk. The limiting distributions 
in this case are given by the generalized central limit theorem , due to the work of Levy, 
Khintchine, Gnedenko and Kolmogorov [17|[27][28 in the 1930. They are called a-stable 
distributions and are presented in appendix [Aj They are characterized by their index of 
stability, a £ (0,2], and include the Gaussian distribution as a special case, for a = 2. 
When a < 2, they have algebraically decaying heavy tails with exponent —(1 + a), and 
infinite variance. For this reason, they are particularly interesting to model stochastic 
processes with high variability, such as solar flare intermittency [29]. By removing the 
hypothesis on the independence of the steps of the CLT, we allow for long-time correlations 
in the Brownian motion. 


In the next section, we introduce the model of the continuous time random walk (CTRW) 
which is a generalization of the random walk. We give two notable examples of this model, 
the Levy-flight and the Levy-walk. Finally, we introduce an other model of non-diffusive 
transport which is based on the Langevin equation of motion and includes long-range 
temporal correlations, the fractional Brownian motion (fBm) and its generalization to a 
non-Gaussian, heavy-tailed process, the fractional Levy motion (fLm). 
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2.1 Continuous time random walk 


The model of the CTRW was first developed in 1965 by Montroll and Weiss to describe 
the mobility of charges in amorphous semiconductors 30 . It has since then found a wide 


range of applications in physics, chemistry, biology, etc. The CTRW supposes that a 
particle, also called random walker, makes successive jumps interrupted by rests. The 
step sizes and the waiting times are drawn from a PDF, called the jump PDF, 

Various choices of the form of il>(x,t ) lead to different situations. For example, if the step 
sizes and the waiting times are independent random variables, the CTRW is decoupled 
and the jump PDF can be written as ip(x,t) = A (x)w(t), where X(x) and w(t) are the step 
size PDF and the waiting time PDF, respectively. In the case of a coupled jump PDF, a 
jump of a certain length involves a certain duration. This is, for example, the case of the 
Levy-walk. 


The CTRW can be described by the equation [311 

/ +oo r+oo 

dx 1 / &t'r)(x\ t / )^(x / — x, t — t') + S(x)S(t), (16) 

-oo J 0 

which links the PDF rj(x,t) of arriving at position x at time t with the PDF ? j(x',t') of 
being arrived at position x' at time i! with a delta Dirac initial condition. 


For the decoupled case, the PDF of the density of walkers is therefore given by 

n(x,t) = f dt'rj(x, t / )T(f — t') dt', (17) 

Jo 

where \k(t) = 1 — Jq w(t') dt' is the probability for a walker of making no jumps between 
the time interval (0, t). 

By taking the space Fourier transform and time Laplace transform^ of eq. ( |l7| ) , one finds 
the Montroll-Weiss equation 


h(k,s ) 


1 — w(s) 
s(l — w(s)X(k )) 


(18) 


Starting from this equation, the fluid limit, meaning that all details of the CTRW that 
are irrelevant at very large temporal and spatial scales are neglected, is usually sought to 
find a equation describing the temporal evolution of the density of random walker. 

For example, classical diffusion is found by taking a Gaussian distribution of step sizes and 
a Poissonian wating time PDF. In fact, the fluid limit is always found to be the classical 

1 In the following we us the notation / for both the Fourier and the Laplace transforms of the function 
/. The difference between the two transforms is indicated by their conjugate variables: x k and 
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diffusion, as long as the characteristic waiting time 


and the variance of step sizes 


r»+oo 


r = 


tw(t) df 


( 19 ) 


A 2 



( 20 ) 


are finite 18 . In this case the diffusion coefficient is given by D = A 2 /r. 


To describe non-diffusive transport, A and w have to be chosen such that r or A 2 , or both, 
diverge. A natural choice is to draw both functions in the family of stable distributions 
(with certain restrictions) since they are the limit of sums of random variables |32|. The 
non-Gaussian property will arise from a divergent A 2 and the non-Markovianity from 
a divergent r. In the case where both characteristic scales diverge, we can choose a 
strictly symmetric stable distribution with characteristic exponent a < 2 for the step sizes 
distribution and a one-sided stable distribution with characteristic exponent j3 < 1 for the 
waiting time distribution (see appendix |A| : 

\(k) = e _0 '“l fe l a ~ 1 — <j a \k\ a for k —»• 0 with a <2 (21) 


and 

w(s) = e~^ sP ~ 1 — for s —> 0 with (3 < 1, 

with the following asymptotic behavior 

A(x) ~ aT (Q+1) , when \x\ —> oo 


( 22 ) 


(23) 


and 

w(t) ~ when t —> oo. 

Using those expression in eq. ( |18| ) and keeping terms up to first order, one finds 

s /3— 1 

n(k,s ) = 




(24) 


(25) 


which can be rearranged as 

s^h(k, s) — s 13 ^ 1 = —°^n\k\ a h(k, s). (26) 

If a = 2 and (3 = 1, one recognizes the Laplace transform of the first order time derivative 
of h(k, t) on the left hand side (lhs) and the Fourier transform of the second order space 
derivative of h(x,s ) on the right hand side (rhs). Therefore, in this case, we recover the 
classical diffusion equation (eq. ([ 7 ])) with D = o 2 1[i. 


In the case were a < 2 and (3 < 1, the lhs and rhs correspond to generalizations of the 
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differential operator to fractional orders (see appendix [5]) 

Q0 9“ 

-^n(x, t) = D a ,p—-n(x, t), (27) 

where D a ^ = o a /pfi { m“/s^] is a generalized diffusion coefficient. This equation is called 
the space and time fractional diffusion equation and describes the time evolution of the 
fluid limit of a CTRW with long-term memory (non-Markovian) and long-range spatial 
correlations (non-Gaussian). When a —> 2 and /5 —>• 1, the classical diffusion equation (eq. 
0 ) is recovered. We call a the spatial transport exponent and (3 the temporal transport 
exponent. The time derivative operator is the Caputo fractional derivative and the space 
derivative is the Riesz fractional derivative (appendix [B]). 

By using the scaling properties of the Fourier and Laplace transform 


/F[f(ax)] = |a| 1 f(k/a) and C[f(bu)] = b l f(s/b), b > 0 , 


-l. 


(28) 


on eq. (25), the following scaling property of the propagator of eq. (27) is inferred 

G at p(x,t) = t~ H K a fi (^ 7 ) , 


(29) 


where H = /3/a is the self-similarity index also called Hurst exponent [33], K a ^ is called 
the reduced Green function and ffr is the similarity variable 20 . This scaling implies that 
the stochastic process associated with the propagator is self-similar with index H = (3/a 
(or iL-self similar). Indeed, if we rescale the time by a factor A > 0, we find 


G a fi{x,\t) = A H t H K a ,g G »A X ^ 

implying that the motion is invariant under the following transformation 


(30) 


(x, t ) i-)- ( X H x , At). 


(31) 


The self-similarity of the process leads then to the following scaling of the moments of the 
distribution with time 

(|x| s )oct^. (32) 

Indeed, if we assume (|x| s ) oc t K , we have 


(lA^xj 8 ) oc (At) 


K 


\ sH (|a;| s ) oc A K t 


K.K 


K = sH. 


(33) 


By identifying eq. (15) with eq. (32) we see that 7 = 2 H. The transport is therefore 


superdiffusive if 2/3 > a and subdiffusive if 2/3 < a. For the special case 2/3 = a, (3 / 
l,a / 2 the transport is called quasidiffusive \ 1 1 and exhibits the same scaling as the 
classical diffusion despite the crucial difference of a non-Gaussian PDF and non-Markovian 
time increments. We would like to point out that only verifying the time dependence of 
the variance of displacements is not sufficient to conclude that a transport process is 
diffusive. Whenever possible, the temporal evolution of higher moments, or of the full 
distribution, should be examined to determine the presence of non-Gaussian features. 
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Figure 2 | Example of a Levy flight trajectory with a = 1.5. Contrary to the Brownian motion 
(%• 0. arbitrary long steps can arise, on all scale length, leading to the clustering nature of the 
Levy flight. 


Evaluating moments of the distribution can be delicate since the integrals that define 
them do not always converge. For example, the s-order moments of a a-stable distributions 
converges only for s < a [l7j. However, in practice, the temporal scaling can be recovered 
by taking “truncated” moments or by computing moments of fractional order smaller than 


a 


18 


2.1.1 Levy flight 


The Levy flight is a particular case of the CTRW with a finite characteristic waiting time, 
t < oo, but step size distribution given by a symmetric stable distribution with diverging 
variance, A 2 —> oo. The trajectories of Levy flights have been shown to model the foraging 
motions of many living organisms [7,34 


Mandelbrot also used this model to simulate the 


fractal galaxy distribution in the Universe 23 . In fact, a fractal dimension D = a can be 


assigned to the trajectories. An example of a Levy flight trajectory is shown in fig. [2] 

The Levy flight can be modeled by taking a Poissonian distribution for the waiting time 
PDF w(t) = T~ 1 e~ t / r , with Laplace transform w(s ) = (1 + sr )^ 1 ~ 1 — ts , for s — > 0. 


Using this in eq. (18), we find 


(jq, (Ay s ) 


s+9\k\ c 


(34) 


which, upon Fourier-Laplace inversion, shows that the propagator of the Levy flight is a 
symmetric stable distribution 


G(x,t) = t H L a , a , 


(35) 


with a self-similar index H = 1/a. Rearranging the terms of eq. (34) and Fourier-Laplace 
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inverting it, we find the space fractional diffusion equation 


describing the Levy flight. 


_9 
9 1 


n(x, t ) = D, 


9 " 

9|x|“ 


n(x,t ) 


(36) 


The Levy flight results in a super diffusive process (with the exception of the case a = 2) 
with a diverging mean-square displacement (x 2 ) —> oo for a < 2. 


The presence of arbitrary long jumps without any restriction on the step duration leads 
to rather unphysical situations 118,19 . One way of solving this is introduced in the Levy 


walk model, which is often more appropriate to describe physical systems. 


2.1.2 Levy walk 


Similarly to the Levy flight, the Levy walk model maintains a diverging variance of dis¬ 
tribution of step sizes. However, a coupling between the step sizes and the step duration 


is included in the jump pdf such that 131 35 


if(x, t) = A(x)<5(|x| — vt v 


(37) 


where A(x) —> \x\ ^ as |x| —> oo and v is a generalized velocity which penalizes long jumps 


such that the variance is finite 19 . Depending on the values of the two exponent p and v 


the transport can be either superdiffusive or subdiffusive. Due to the coupled form of the 
jump pdf, the derivation of a transport equation describing the evolution of the PDF has 
only been achieved recently in the case v = 1, by using a fractional version of the material 


derivative 36 . 


2.2 Langevin approach 

A microscopic description of Brownian motion equivalent to the one of Einstein presented 
in sec. [I] was introduced by Langevin in 1908 116] An uncorrelated Gaussian noise, rep¬ 
resenting the random force due to the interaction with the fluid molecules, is used in the 
equation of motion of a test particle. The equation of motion becomes a stochastic equa¬ 
tion, whose average motion shows the same diffusive scaling. Fractional Brownian motion 
(fBm) introduces long-range temporal dependence in the Gaussian noise, which can lead to 
a non-linear scaling of the positional variance. On the other hand, non-Gaussian statistics 
can be introduced by choosing a non-Gaussian noise. For example, stable Levy motion |l7| 
replaces the Gaussian noise with a random noise distributed according to a Levy stable 
distribution with heavy-tails (appendix [A|. 

The classical Langevin equation is written [16 

mx(t) = ~m'yx(t) + f(t), (38) 
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where m is the mass of the test particle, 7 is the friction coefficient and £(f) is the random 
force due to the random collisions with the surrounding particles. In the case of the 
Brownian motion, £(f) is a a white noise, i.e. a Gaussian noise with an infinitely short 
correlation time: (£(fi)£(f 2 )) = C8{t\ — £ 2 ), where C is a constant. In eq. 38. the forces 
acting on the particle are separated in two groups, the macroscopic, slowly varying ones, 
represented by the dissipative force —m'yv, and the microscopic, rapidly varying ones, 
represented by the fluctuating force £(f). 

If the time scale of the particle motion is comparable to the time scale of the collisions, 
the assumption of a white noise and a constant friction have to be abandoned. This leads 
to the generalized Langevin equation (GLE) |37] 


x(t) = — / /3(f — df r + £(f), 

Jo 


(39) 


where m = 1 is used for simplicity. Here, /3(f) is the memory kernel and £(f) is the random 
force which is zero-centered and stationary, i.e. (£(fi)£(f 2 )) = C(|fi — ^ 2 1) = C'(r). The 
fluctuation-dissipation theorem 37-40 states that the dissipation is the macroscopic man¬ 
ifestation of the disordering effect of the fluctuations and relates the correlation function 
of the random forces C(t ) with /3(f) by 


k B TP(t) = C(t). 


(40) 


Assuming x(0) = 0, u(0) = i>o and Laplace transforming eq. (39), we find 

vo + (3(s) +i(s) 


x(s) = 


s(s + (3(s)) 


(41) 


where /3(s) and £(s) are the Laplace transforms of /3(f) and £(f). Upon Laplace inversion, 
one finds the equation of the particle position 


x(t) = v 0 H(t)+ / H{t - df', 


(42) 


where H(t) is the relaxation function 40 defined by its Laplace transform 

H{s) = 1 


s(s + /3(s)) 


(43) 


We note that, in accordance with the classical Langevin equation, if /3(f) = 7 = cste., 
the relaxation function, H{t) = -(1 — e -7i ), is exponentially decaying and the position is 
given by 


x(t) = 52(1 - e -7 + - f 
7 7 ./o 


(1 


_ e -7(t-t'))^(f')dt / . 


For f —> 0 the ballistic motion x(t) = vot is recovered and for f ^ 7 1 , 

nt 


x(t) = — + - f £(t')dt'. 
7 7 Jo 


(44) 


(45) 
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When vq = 0, eq. (45) is referred to as ordinary Brownian motion (oBm). 


If we take £ as a white noise and (vq) = 0, the variance of the displacement is given by 

{x 2 {t)) = ^!^J €{t-h)dtij £(t-t 2 )dt 2 ^ 

= £(/ dh J dt 2 £(t - ti)£{t - t 2 )^ 

= \ { dti [ dt 2 (£(t - ti)£(t - t 2 )) 

= 4 [ dti [ dt 2 Cd(t 2 — ti) 

T Jo Jo 

= 4 [ t Cdt 1 = ^t= —t, (46) 

T Jo T 7 


where we used eq. (40). Thus, we recover the linear temporal scaling of the mean-squared 
displacement of the classical diffusion (eq. ©)■ 


2.2.1 Fractional Brownian motion 


Fractional Brownian motion (fBm) was proposed by Mandelbrot and Van Ness in 1968 |41 
to model the variations of cumulated water flows in the great lakes of the Nile river basin 
observed by Hurst 1421. Hurst studied the record of river level and other physical quan¬ 
tities such as rainfall, temperature, pressure, the growth of tree rings, sunspot numbers 
and wheat prices. He found that the range of those records, rescaled by their standard 
deviation, is proportional to t H , where t is the time and 1/2 < H < 1 is, ever since, called 
the Hurst exponent. Since then, it has found a wide range of applications in systems 
showing long time interdependence. 


Slightly different representations exists in the literature, here we use the following 143 


x H (t) = 


T(H + 1/2) J o 


(t - t') H ~ 1/2 t(t') dt', 


(47) 


where xn(t) represents the position of a particle experiencing fBm, £(i) is a Gaussian 
uncorrelated noise, H £ (0,1] and T(-) is the gamma function. FBm is constructed 


as a moving averaged of the ordinary Brownian motion (oBm) (eq. (45)), in which past 
increments are weighted by the power law kernel . It has a zero mean (x#(t)) = 

0 for H < 1. From its definition and the fact that the Gaussian noise is self-similar with 
exponent 1/2, it follows that is H -self similar (eq. ©) and that it has stationary 

increments, xn(t) — xjj(s) = xn(t — s) [17]. Using these two properties, we can show that 
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the correlation function is 


(x H (t)x H (s)) = K h ^ {(x 2 H {t)) + (x 2 h (s)) - <( x H (t ) - ®^(s)) 2 )} 

= -Ktf ^ {(4W) + ( x h{s)) ~ (:x 2 H {t - s))} 

= K H l -{t 2H + s 2H -\t-s\ 2H ) , (48) 

where I\h is a positive constant and where we recall that for a iL-self similar process, the 
variance scales with time as 

<*h(0> « t2H ■ ( 49 ) 

We note that fBm is subdiffusive for 0 < H < 1/2, superdiffusive for 1/2 < H < 1, 
ballistic for H = 1 and correspond to oBm for H = 1/2. 


The increments of the fBm, £#(t), is a stationary Gaussian process known as fractional 
Gaussian noise (fGn) and defined as 


x H (t) 


6/(0 d t'. 


(50) 


The correlation function, Cn(t), of £#(t) is given by the derivative of equation (48) with 
respect to t and s 


C H (\t-s\) = (&K*)6r00> = 2K H H{2H-l)\t-s\ 2H - 2 + 2K H H\t-s\ 2H - 1 5(t-s). (51) 


We note that Cn(\t — s|) behaves as a power law for t = \t — s| — > oo and recovers the 
ordinary Brownian behavior, Cn(\t — s|) = Ki/ 2 6(t — s), for H = 1 / 2 . 

The function Ch{t) tends to zero for r = \t — s| — > oo for 0 < H < 1, but when 
1/2 < H < 1 , £#(i) exhibits long-range dependence, i.e. Ch(t) tends to zero so slowly 
that /I/ 0 Ch(t)cIt = oo. It is said to be correlated. For 0 < H < 1/2, there is no long-range 
dependence, but the coefficient (2 H — 1) is negative [lT . In this case the £/j(t) is said to 
be anti-correlated. Figure [3] shows three examples of fBm trajectories. 


In the framework of the generalized Langevin equation (GLE), it is possible to find the fBm 
as a solution by using a random force with long-range correlations, namely with a power- 
law correlation function. The memory kernel, /3(t), is then found with the fluctuation- 


dissipation theorem (eq. (40)) and, consequently, also have a power-law form. When the 
random force is chosen to be the fractional Gaussian noise (fGn), the GLE can be written 
as a fractional differential equation (44|, however, the solution of this equation is limited 
to the subdiffusive and diffusive case. From the physical point of view, the superdiffusive 
case is found only when the random force is “external”, meaning that the fluctuation- 
dissipation theorem does not hold and that the driving noise and the dissipation may have 
different origins, which may be the case in nonequilibrium systems 40 . 


It is possible to find the propagator of the fBm by using the method of path integrals 43 


borrowed from quantum mechanics (note that Shrodinger’s equation resembles a diffusion 
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Figure 3 | Examples of realizations of the fractional Brownian motion (fBm) for different values 
of the self-similar index H. For H = 0.8, the increments of the trajectory are positively correlated 
which results in a persistent motion. The oBm with uncorrelated increments is retrieved for 
H = 0.5 and for H = 0.2, the increments of the trajectory are anti-correlated which results in a 
anti-persistent motion. 


equation with an imaginary diffusion coefficient). The propagator, given by 


^ ^ [HT(H +1/2) 

G H (x,t) = \j — -^-exp 


-HT 2 {H+ l/2) w 


(52) 


is a Gaussian function with a variance proportional to t 2H . It has the same form than 
the propagator of the oBm (eq. (|l3|)) but with a “stretched” time t 2H . The transport 


equation of the fBm is easily found from the Fourier transform of eq. (52) to be 


0 0 ^ 

—n(x,t) = Dpt 13 - 1 —n(x,t), 


(53) 


45 


• By 


where 0 < /3 = 2 H < 2 and Dp = (2 T 2 (H + 1/2)) -1 is a stretched diffusion coefficient of 
dimensions [m 2 /s l3 ]. This equation is called the stretched time diffusion equation 
using the rule 

a „ . a 

(54) 


3th 9t’ 


it can be interpreted as the result of the classical diffusion equation with a stretched time. 
We note that the Langevin approach and the CTRW approach are not equivalent in the 


non-Markovian case. Equation (53) is local in time, whereas eq. (27) with a = 2, /3 < 1, the 


time-fractional diffusion equation , is not. In the fBm case, the non-Markovian character 
is provided by a time dependent diffusivity D = Moreover, The solution of eq. 


(53) is Gaussian, while the solution of the time-fractional diffusion equation is not; it is 


given by the transcendental functions known as the M-Wright function which tends to the 
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Gaussian function for f3 = 1 


45 


2.2.2 Fractional Levy motion 


Here, we discuss the fractional Levy motion (fLm), which is a generalization of the fBm, 
including both long-range temporal dependence and non-Gaussian statistics. 


The stochastic equation defining the fLm process is |46 47 

rt 


X a ,H{t ) = 


1 


T (H — l/a + l) 


'o 




(55) 


where £o, i(r (t) is an uncorrelated noise distributed according to a Levy symmetric, strictly 
stable distribution, with index of stability a, (0 < a < 2) and scale parameter a. From 
the properties of a-stable random variable, we have x a = A H x a ^n(t), with H = 

/3/a. Therefore, the fLm belongs to the important family of H -self similar process with 
stationary increments (also abbreviated H-sssi ), like the fBm. Consequently, the moments 
of Xa^nit) exhibit the desired general non-classical feature 


( \x a ,H{t )\ s ) oc t sIi , 


(56) 


where 0 < s < a, to ensure convergence of the moments. For a non-degenerated process, 
the values of H are restricted to [l7j 


0 < H < 1/a if a < 1, 
0 < H < 1 if a > 1. 


(57) 


The long-term memory is engendered by the convolution with the power-law kernel and the 
non-Gaussian statistics by the Levy noise. The fLm generalizes the fractional Brownian 
motion (fBm) |41 . Indeed, for a = 2, the noise has a Gaussian distribution and the process 
is the fBm. When H = 1/a the process is time-uncorrelated and when H < 1/a or H > 
1/a the process exhibits negative or positive time correlations, respectively. Therefore, for 
a = 2 and H = 1/2, one recovers the oBm corresponding to classical diffusion (eq. ©)• 


Using path integrals, Calvo, Sanchez and Carreras have shown that the transport equation 
of the fLm process is a space-fractional diffusion equation with time dependent diffusivity 

0 

r\ 

ai n(x,t) = Kt 13 - 1 —-^n(x,t). (58) 

ot d \x\ 

Here n(x, t ) is the density of particles, K is an effective diffusion coefficient and a and 
/3 are the space and time transport exponents, respectively, with H = (3/a. The space 
derivative of order a is the Riesz fractional differential operator [48] (appendix [B]). The 
restriction on the range of permissible values for H (eq. ©) translates for (3 as 


0 < /3 < 1 if a < 1, 

0 < /3 < a if a>l. 


(59) 
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Figure 4 | Examples of propagators for different parameters of the fLm for times t-i = 1,2, ...,9 
and effective diffusivity K = 1. a: Subdiffusive case with a = 2 (Gaussian) and /3 = 0.5 (anti¬ 
correlated) (H = 0.25). b: Diffusive case with a = 2 (Gaussian) and /? = 1 (uncorrelated) 
( H = 0.5). c: Superdiffusive case with a = 1 (Levy stable) and /3 = 1 (uncorrelated) (H = 1). 


When /3 < 1 or /? > 1 the process is negatively or positively time correlated, (3 = 1 
corresponding to an uncorrelated process. 


The propagator of eq. (58) is a Levy distribution which depends on x/t 


H 


47 


o *\ _ C ( a >P) T 

~ +f}/a 


C(a,P) 


x 


tP/ a \ ’ 


(60) 


where C(a, /3) = /3‘/« r + l) and a = K 1 / a T + l). 


Again, for a = 2 the space-fractional derivative becomes a second order derivative and the 


propagator is a Gaussian corresponding to the case of the fBm (eq. (52) and eq. (53)). 


For a = 2 and (3 = 1 (H = 1/2), eq. (58) becomes the classical diffusion equation and 


the propagator has the well known form of a Gaussian with variance growing linearly with 
time (eq. ([7]) and eq. ©). Finally, for [3 = 1 and a < 2, eq. ( |60| and eq. ( [58] ) are 
Markovian and correspond to the propagator and transport equation of the Levy flight 
(eq. (35) and eq. (36)). Examples of fLm propagators are shown in eq. 0 and the 
different transport regimes of the fLm, as a function of a and (3 are summarized in fig. [5j 

We note here that a generalization of the model of the fLm to include asymmetric proba¬ 
bility density functions as recently been derived and applied to the case of suprathermal 
ion transport in turbulent plasmas |9,13 . 
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Figure 5 | Parameter space for the fLm. The values of H = 0/a lie in the shaded region and 
depending on the values of the spatial transport exponent, a, and the temporal transport exponent, 
0, the transport can be superdiffusive (H > 1/2) or subdiffusive ( H < 1/2). Diffusive transport is 
found for H = 1/2 (white line), the classical diffusion corresponding to the case 0 = 1 and a = 2 
(black dot). Gaussian transport is found for a = 2 and Markovian transport for 0 = 1. For 0 > 1 
the motion is persistent and for 0 < 1 it is anti-persistent. The fBm (green line) is found for a = 2 
and 0 < 0 < 2 and the Levy flight (yellow line) for a < 2 and 0 = 1. A similar figure can be 
drawn for the CTRW, but with the difference that the non-Markovianity is not due to persistence 
or anti-persistence but to the non-locality in time and, in this case, 0 < 1. 


A Stable distributions 


We recall here the main properties of a-stable random variables and distributions (also 
called Levy a-stable distributions), although we refer the reader to reference [IT] for a 
detailed monograph. The generalized central limit theorem states that stable distributions 
are the limiting distributions of normalized sums of independent, identically distributed 
random variables. The Gaussian distribution is a particular case that corresponds to the 
case where the random variables have a finite variance. Stable distributions are especially 
interesting because, with the exception of the Gaussian distribution, they allow heavy- 
tails and non-zero skewness. Their probability densities exist and are continuous, but 
their expressions in closed form exist only for a few particular cases (Gaussian, Cauchy 
and Levy distributions) and they are usually described by their characteristic function, 
i.e. their Fourier transform. 

A random variable X is said to be stable if for any A, B > 0 there is C > 0 and l)£l 
such that [TF| 

AX + BX = CX + D. (61) 

Moreover, for any stable random variable there is a number a 6 (0, 2] such that 

C a = A a + B a . (62) 
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Figure 6 | Example of stable distributions on a lin-lin (left) and log-log scale (right). All the 
distributions have shift parameter p = 0 and scale parameter a = 10. A Gaussian distribution is 
shown in black (a = 2), a distribution with a = 1.2 in blue and a one sided distribution (totally 
skewed to the right) in red. The log-log scale shows the heavy tails of the distributions with a < 2, 
with exponent — (a + 1). 


The Gaussian distribution corresponds to the case a = 2. 

Different parameterizations of their characteristic function are possible. We adopt here 
the following parameterization for the stable distribution L e a a (close to parameterization 
C in (49)): 

^[ Ld a,aAx)}(k) = exp |—o- Q |fc| Q e® slgn(fc) ^ + ikp} . (63) 

Here a E (0,2] is the index of stability or characteristic exponent , a > 0 is the scale 
parameter, 8 is the asymmetry parameter (|0| < min(a,2 — a)) and p E M the shift 
parameter. For example, the Gaussian distribution is found for a = 2 and has mean 
equal to p and standard deviation equal to y/2a (the asymmetry parameter 8 is equal to 
zero when a = 2). When p = 0, the class of distributions reduces to the strictly stable 
distribution L e a a and when p = 8 = 0, the class of distributions reduces to the symmetric 
a-stable distributions, L a ^ a . 

When a < 2, stable distributions have the interesting property of having one tail, if their 
skewness is maximum, or both tails that behave asymptotically as power laws (heavy- 
tail) 117], 

L a,*A x ) ~ , when |x| -A oo. (64) 

As a result, they always have infinite variance when a < 2, which reflects their capability 
of modeling processes with large fluctuations. For a < 1 they also have infinite first 
moments. 

Stable distributions that have a < 1 and are totally skewed to the right or to the left 
( 8 = ±a) are one sided. They are only defined for x > 0 if 8 = —a and for x < 0 if 8 = a. 
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For /j, = 0 and 6 = —a, their Laplace transform is given by jTr] 

£-{L a ,cr{x)\{s) = exp /-(65) 

l cos 2 j 

Although no expression in closed form exists for the general stable distributions (they are 
expressed with transcendental functions known as Mittag-Leffler functions (201), numerical 
algorithms allow to compute their probability distribution function with great accuracy 
and very efficiently (see [501 for example). 


B Fractional differential operators 


The idea of generalizing the differential operation to fractional order is as old as differen¬ 
tial calculus. Leibniz, Euler, Liouville, Riemann and Fourier are among the many great 
mathematicians who developed fractional differential calculus. It is only more recently, in 
the second half of the 20th century, that it began to be applied in physics and engineering, 
to problems such as the modeling of viscoelasticity in materials 1511 or dynamical processes 
in fractals 1521. 


Among the different definitions, one of the most famous is the Riemann-Liouville differ¬ 


ential operators that can be defined explicitly by means of the integral operators 32,53 
The left and right Riemann-Liouville fraction derivative of order a are 


D £/(*) = 


1 


d n 


T(m — a) dx^ 


f(x') 


(;x - x ') a - m+l 


dx 


( 66 ) 


and 

b D £/(*) 

where T(t) = j^ ’ x* 1 e x dx is the gamma function, m is the integer satisfying m — 1 < 
a < m. and a and b are the start and end point of the operators. We immediately see from 
there definitions that an important difference of the fractional version of the derivative of 
a function at a point x is that it is not a local property. As a matter of fact they can 
depend on the value of the function very far from x. In the cases in which the start point 
a or the end point b extend all the way to infinity, the following notation is generally used 


-1 


r(m — a) d(— x) r 


f{x') 


_ 1 


(x — x') 


dx 


(67) 


d a f 

dx" 


= D" 

-oo *-> x 


fix) 


and 


d °f 
d(—x)" 


_~hoo 


D “/(*). 


( 68 ) 


The Fourier transform of these operators sheds light on their signification and their usage as 
they appear as a natural generalization of the Fourier transform of the derivative operator 


T 


' d a r 

dx" 


(■ ~ik) a f(k ) 


and 


T 


d a f 

d(—x)" 


m a f(k). 


(69) 


A symmetrization of these operators leads to the Riesz fractional derivative operator 
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-1 


(70) 


d a f 

d|x|“ 


2 cos(7ra/2) 


with the following Fourier transform 


d°7 d a f \ 
Ax a d(— x) a ) 


T 


d Q / 
d|a;| a 


ifcr/>). 


(71) 


A different definition than the Riemann-Liouville is the Caputo fractional derivative op¬ 


erator of order f3 20,51 


d p t m = 


r(m — (5) 




I o (t-t')P~ m+1 


dt 


with m — 1 < (3 < m. 


(72) 


This definition of the fractional derivative is usually associated with derivatives in time 
because of the practical form of its Laplace transform 


m -1 jfc f 

= < 73 ) 

k =0 

which depends only on the initial values of f(t) and its integer derivatives. The Laplace 
transform of the Riemann-Liouville derivative depends instead on the initial values of 
the fractional derivatives of lower order than f3, which makes is it not practical for real 
applications 1511. 


£ 
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